{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "32fa5833",
   "metadata": {},
   "source": [
    "# Magnetostatics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fdb3307",
   "metadata": {},
   "outputs": [],
   "source": [
    "from netgen.occ import *\n",
    "from ngsolve import *\n",
    "from ngsolve.webgui import Draw\n",
    "from netgen.webgui import Draw as DrawGeo\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc952d2a",
   "metadata": {},
   "source": [
    "### model of the coil:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10a2b96b",
   "metadata": {},
   "outputs": [],
   "source": [
    "cyl = Cylinder((0,0,0), Z, r=0.01, h=0.03).faces[0]\n",
    "heli = Edge(Segment((0,0), (12*math.pi, 0.03)), cyl)\n",
    "ps = heli.start\n",
    "vs = heli.start_tangent\n",
    "pe = heli.end\n",
    "ve = heli.end_tangent\n",
    "\n",
    "e1 = Segment((0,0,-0.03), (0,0,-0.01))\n",
    "c1 = BezierCurve( [(0,0,-0.01), (0,0,0), ps-vs, ps])\n",
    "e2 = Segment((0,0,0.04), (0,0,0.06))\n",
    "c2 = BezierCurve( [pe, pe+ve, (0,0,0.03), (0,0,0.04)])\n",
    "spiral = Wire([e1, c1, heli, c2, e2])\n",
    "circ = Face(Wire([Circle((0,0,-0.03), Z, 0.001)]))\n",
    "coil = Pipe(spiral, circ)\n",
    "\n",
    "coil.faces.maxh=0.2\n",
    "coil.faces.name=\"coilbnd\"\n",
    "coil.faces.Max(Z).name=\"in\"\n",
    "coil.faces.Min(Z).name=\"out\"\n",
    "coil.mat(\"coil\")\n",
    "crosssection = coil.faces.Max(Z).mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01d1ff33",
   "metadata": {},
   "outputs": [],
   "source": [
    "DrawGeo (coil);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed4c0f9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "box = Box((-0.04,-0.04,-0.03), (0.04,0.04,0.06))\n",
    "box.faces.name = \"outer\"\n",
    "air = box-coil\n",
    "air.mat(\"air\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73262a2c",
   "metadata": {},
   "source": [
    "### mesh-generation of coil and air-box:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7a5847a",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "geo = OCCGeometry(Glue([coil,air]))\n",
    "with TaskManager():\n",
    "    mesh = Mesh(geo.GenerateMesh(meshsize.coarse, maxh=0.01)).Curve(3)\n",
    "    \n",
    "Draw (mesh, clipping={\"y\":1, \"z\":0, \"dist\":0.012});"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef85e116",
   "metadata": {},
   "source": [
    "checking mesh data materials and boundaries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d349c6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh.ne, mesh.nv, mesh.GetMaterials(), mesh.GetBoundaries()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c9b3192",
   "metadata": {},
   "source": [
    "### Solve a potential problem to determine current density in wire:\n",
    "on the domain $\\Omega_{\\text{coil}}$:\n",
    "\\begin{eqnarray*}\n",
    "j & = & \\sigma \\nabla \\Phi \\\\\n",
    "\\operatorname{div} j & = & 0\n",
    "\\end{eqnarray*}\n",
    "port boundary conditions: \n",
    "\\begin{eqnarray*}\n",
    "\\Phi & = & 0  \\qquad \\qquad \\text{on } \\Gamma_{\\text{out}},  \\\\\n",
    "j_n & = & \\frac{1}{|S|} \\quad \\qquad \\text{on } \\Gamma_{\\text{in}},\n",
    "\\end{eqnarray*}\n",
    "and $j_n=0$ else"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce0f3fe4",
   "metadata": {},
   "outputs": [],
   "source": [
    "fespot = H1(mesh, order=3, definedon=\"coil\", dirichlet=\"out\")\n",
    "phi,psi = fespot.TnT()\n",
    "sigma = 58.7e6\n",
    "bfa = BilinearForm(sigma*grad(phi)*grad(psi)*dx).Assemble()\n",
    "inv = bfa.mat.Inverse(freedofs=fespot.FreeDofs(), inverse=\"sparsecholesky\")\n",
    "lff = LinearForm(1/crosssection*psi*ds(\"in\")).Assemble()\n",
    "gfphi = GridFunction(fespot)\n",
    "gfphi.vec.data = inv * lff.vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a353b78",
   "metadata": {},
   "outputs": [],
   "source": [
    "Draw (gfphi, draw_vol=False, clipping={\"y\":1, \"z\":0, \"dist\":0.012});"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "453d1083",
   "metadata": {},
   "source": [
    "### Solve magnetostatic problem:\n",
    "\n",
    "current source is current from potential equation:\n",
    "\n",
    "$$\n",
    "\\int \\mu^{-1} \\operatorname{curl} u \\cdot \\operatorname{curl} v \\, dx =\n",
    "\\int j \\cdot v \\, dx\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9fe77bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = HCurl(mesh, order=2, nograds=True)\n",
    "print (\"HCurl dofs:\", fes.ndof)\n",
    "u,v = fes.TnT()\n",
    "mu = 4*math.pi*1e-7\n",
    "a = BilinearForm(1/mu*curl(u)*curl(v)*dx+1e-6/mu*u*v*dx)\n",
    "pre = Preconditioner(a, \"bddc\")\n",
    "f = LinearForm(sigma*grad(gfphi)*v*dx(\"coil\"))\n",
    "with TaskManager():\n",
    "    a.Assemble()\n",
    "    f.Assemble()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a950e839",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ngsolve.krylovspace import CGSolver\n",
    "inv = CGSolver(a.mat, pre, printrates=True)\n",
    "gfu = GridFunction(fes)\n",
    "with TaskManager():\n",
    "    gfu.vec.data = inv * f.vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "616c675e",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "Draw (curl(gfu), mesh, draw_surf=False, \\\n",
    "      min=0, max=3e-4, clipping = { \"y\":1, \"z\" : 0, \"function\":False}, vectors = { \"grid_size\":100});"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b7491f4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
